;+
; NAME:
;       OMNI_MASS_COMPLETENESS
;
; PURPOSE:
;       Computes the mass completeness function for use in correcting
;       the mass function at the low end.
;
; CATEGORY:
;       distance-omnibus Mass Routine
;
; CALLING SEQUENCE:
;       complete_fn = OMNI_MASS_COMPLETENESS (mscale,  [SAMPLETD=sampletd])
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       CONFFILE -- Name of the configuration file to use for survey
;                   information [Default: conffiles/survey_info.conf]
;       NREAL    -- Number of realizations of the mass function to
;                   compute [Default: 1000].
;       SAMPLETD -- Override of sampletd value in mass_deriv.conf, may
;                   be 0 or 1.  [Default: value in mass_deriv.conf]
;       DMFBINW  -- Override of dmfbinw value in mass_deriv.conf, in
;                   units of Log10 M.  [Default: value in
;                   mass_deriv.conf]
;
; KEYWORD PARAMETERS:
;       ML      -- Compute the Maximum-Likelihood mass function only.
;                  If set, NREAL is ignored.
;       MEGA    -- Compute a MEGA mass function, with NREAL number of
;                  stacked realizations.  Exclusive of ML -- If both
;                  set, then ML will be chosen.
;       PLOT    -- Create a plot of the mass completeness function.
;       EBARS   -- Compute the error bars on the CFN
;       VERBOSE -- Cause the routine to contract logorrhoea.
;       SILENT  -- Mum's the word.
;
; OUTPUTS:
;       COMPLETE_FN -- The mass completeness function, interpolated to
;                      a regular grid based on DMFBINW.
;       MSCALE      -- The mass scale accompanying COMPLETE_FN.
;
; OPTIONAL OUTPUTS:
;       CLIM90 -- The 90% completeness level.
;       CLIM50 -- The 50% completeness level.
;       CLIM10 -- The 10% completeness level.
;
; COMMON BLOCKS:
;       OMNI_CONFIG -- The set of configuration structures, read in
;                      from the config files in conffiles/
;
; MODIFICATION HISTORY:
;
;       Created:  06/04/13, TPEB -- Initial version.
;       Modified: 06/10/13, TPEB -- Name change to OMNI_* and
;                                   conversion into a function with
;                                   output completeness function on a
;                                   regular grid.
;       Modified: 06/13/13, TPEB -- Hacked out obsolete code between
;                                   GOTO staemements.
;       Modified: 07/17/13, TPEB -- Added 68.3% confidence region in
;                                   completeness function for the MEGA
;                                   cases, to ascertain uncertainty in
;                                   the, say, 90% completeness level.
;                                   Also, some minor code cleanup.
;       Modified: 09/04/13, TPEB -- Added the 10% completeness
;                                   output.
;
;-

FUNCTION OMNI_MASS_COMPLETENESS, mscale, CONFFILE=cfile, ML=ml, MEGA=mega, $
                                 NREAL=nreal, SAMPLETD=sampletd, $
                                 DMFBINW=dmfbinw, CLIM90=clim90, $
                                 CLIM50=clim50, CLIM10=clim10, $
                                 PLOT=plot, VERBOSE=verbose, $
                                 SILENT=silent, EBARS=ebars,$
                                 MASSES=cmfn
  
  COMPILE_OPT IDL2, LOGICAL_PREDICATE
  
  COMMON OMNI_CONFIG, conf, mw, local, dpdfs, ancil, fmt, conffile
  
  ;; Parse keywords
  mega   = KEYWORD_SET(mega) && ~KEYWORD_SET(ml) ; Compute a MEGA mfn?
  ml     = ~mega                ; Compute the Maximum-Likelihood mfn?
  nreal  = ~n_elements(nreal) ? 1000L : long(nreal)
  plot   = KEYWORD_SET(plot)
  ebars  = KEYWORD_SET(ebars)
  ndloop = ml ? 1L : long(nreal) ; used for plotting histograms
  
  ;; Read in survey-info, galactic-params, & dpdf-params config files
  conf     = omni_load_conf(cfile)
  dpdfs    = omni_read_conffile('./conffiles/dpdf_params.conf')
  mconf    = omni_read_conffile('./conffiles/mass_deriv.conf')
  d        = dindgen(dpdfs.nbins)*dpdfs.binsize + dpdfs.binstart
  sampletd = n_elements(sampletd) ? sampletd GE 1 : mconf.sampletd
  
  ;; Read in the SURVEY catalog 
  restore,'./local/'+conf.survey+'_map_locations.sav',VERBOSE=verbose
  
  
  ;; Some physics for a little later...
  h     = 6.626068d-34          ; Planck constant    [m^2 kg s^-1]
  kb    = 1.3806503d-23         ; Boltzmann constant [m^2 kg s^-2 K^-1]
  hnuk  = h*conf.nu/kb          ; h*nu/kb            [K]
  
  ;;=================================================================
  ;; Get the desired mass function by calling OMNI_MASS_FUNCTION with
  ;; the /MFNSAVE option to save the mass function to disk for us to
  ;; read in here.
  OMNI_MASS_FUNCTION, CONFFILE=cfile, NREAL=nreal, ML=ml, MEGA=mega, $
                      SAMPLETD=sampletd, /MFNSAVE, MFNFN=mfnfn, $
                      VERBOSE=verbose, SILENT=silent
  IF verbose THEN print,'IDL save file: ',mfnfn
  restore,mfnfn,VERBOSE=verbose
  
  ;; Limit SURVEY structure to the items saved in the above line
  surv_ind = WHERE_ARRAY(s.cnum, survey.cnum, n_survey)
  IF n_survey NE n_elements(s) THEN message,'Rum now, Think later!!!'
  survey = survey[surv_ind]
  nkp = n_survey
  
  ;; The structure element SURVEY.NOISE contains the mean of the noise
  ;;   map within the source contour.  Use survey.complete x this
  ;;   value as the completeness flux FOR THIS SOURCE.
  smin = survey.noise * conf.complete
  
  ;; For stats, compute the signal:completeness ratio
  scr = s.(iflux) / smin
  
  ;; Plot the histogram of the and completeness fluxes SCR for analysis
  IF plot THEN BEGIN
     myps,'./masses/plots/source_flux_complete.eps',xsize=4.2,ysize=2.1,$
          thick=1.5
     !x.omargin = [-4,-5]
     !y.omargin = [0,-1.25]
     multiplot,[2,1],xgap=0.042,/doyaxis
     csf = 0.55
     symsize = 0.2
     bin = 0.1
     
     ;; Panel 1
     plothist,alog10(smin),bin=bin,charsize=1.0*csf,/noplot,xarr,yarr
     plothist,alog10(s.(iflux)),bin=bin,charsize=1.0*csf,xst=5,$
              color='magenta',/fill,fcolor='thistle',yr=[0,500],$
              ytit='N per '+string(bin,format="(F0.1)")+'-dex bin',$
              xtickformat='exponent10',xr=[-1.5,2.5],yminor=5
     plothist,alog10(smin),bin=bin,charsize=1.0*csf,/over
     cgAxis,xaxis=0,xtit='Flux Density  [Jy]',xr=10.^(!x.crange),$
            /xst,/save,charsize=1.0*csf,/xlog,xtickformat='exponent10'
     cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
     al_legend,/top,/right,linestyle=0,box=0,color=['magenta','black'],$
               ['S!d1.1!n','S!dcomplete!n'],linsize=1.5,charsize=0.9*csf

     multiplot,/doyaxis,/doxaxis
     ;; Panel 2

     xr = [-1.2,0.3]
     cgPlot,smin,s.(iflux),psym=16,symsize=0.4,xtit='S!dcomplete!n  [Jy]',$
            ytit='S!d1.1!n  [Jy]',/xlog,/ylog,charsize=1.0*csf,$
            /nodata,ytickformat='exponent10',xtickformat='exponent10',$
            xr=10.^xr,yr=10.^[-1.5,2.5];,xst=5
     one2one,linestyle=3,color='cyan',/log
     cgOplot,smin,s.(iflux),psym=16,symsize=symsize
     ind = where(scr LE 1.,nn)
     ffrac = double(nn) / n_elements(scr)
     print,'FRAC w/ SCR < 1: ',ffrac
     al_legend,/bottom,/right,charsize=0.9*csf,box=0,$
               ['f < 1:1 = '+string(ffrac,format="(F0.3)")]
     
     ;; ticks = LogLevels(10.^xr)
     ;; ntick = n_elements(ticks)
     
     ;; cgAxis,xaxis=0,xr=10.^xr,/xst,charsize=1.0*csf,/save,/xlog,$
     ;;        xtit='S!dcomplete!n  [Jy]',xticks=ntick-1,xtickv=ticks;; ,$
     ;;        ;; xtickformat='exponent'
     ;; cgAxis,xaxis=0,/xst,xtickformat='blank_axis'
     ;; cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
     
     myps,/done,/mp
  ENDIF
  
  
  
  ;;============================================================
  ;; Compute the mass completeness function, and make plots, as
  ;;   specified.
  
  ;; Place information from OMNI_MASS_FUNCTION and SURVEY into more
  ;;   useful forms for our purposes.
  objs.dist /= 1.d3
  dind = where(objs.dist GE mconf.dmin AND objs.dist LE mconf.dmax, nd)
  darr = (objs.dist)[dind]
  sarr = (cmreplicate(smin,ndloop))[dind]
  
  tarr = mega ? (objs.temp)[dind] : mw.td
  tempcorr = (exp(hnuk/tarr)-1.d) / (exp(hnuk/mw.td)-1.d)
  
  ;; Compute the 'completeness mass function' using the
  ;;   'minimum complete flux' for each source (smin)
  cmfn = m0 * sarr * darr * darr * tempcorr
  cdf = ccdf(cmfn,/cdf)         ; Cumulative distribution function
  
  ;; For the output function, use points at 0.1 x binw, and span from 2
  ;;   orders of magnitude below min(mfn) to 2 orders of magnitude
  ;;   above max(mfn)
  binw = mconf.dmfbinw NE 0. ? mconf.dmfbinw : fdr_width(alog10(cmfn))
  mincfn = floor( alog10(min(cmfn)) - 2. )
  maxcfn = ceil(  alog10(max(cmfn)) + 2. )
  nptcfn = long((maxcfn - mincfn) / (binw / 10.) + 1)
  
  mscale = 10.^(findgen(nptcfn) * (binw / 10.) + mincfn) ; Mass scale
  nbin   = n_elements(mscale)
  
  ;; Compute the completeness function, limited by [0,1] to deal with
  ;;   extrapolation effects in interpol(), plus completeness levels.
  cfn = (interpol(cdf,cmfn,mscale) > 0.) < 1.
  clim90 = interpol(mscale,cfn,0.9) ; The 90% completeness level
  clim50 = interpol(mscale,cfn,0.5) ; The 50% completeness level
  clim10 = interpol(mscale,cfn,0.1) ; The 10% completeness level
  
  ;; If desired, compute the error regions on the cfn for MEGA
  IF mega && ebars THEN BEGIN
     
     e_cfn = fltarr(nreal,nbin)
     ;; Go through each realization
     FOR ii=0, nreal-1 DO BEGIN
        d        = objs[ii].dist
        dind     = where(d GE mconf.dmin AND d LE mconf.dmax, nd)
        darr     = d[dind]
        sarr     = smin[dind]
        tarr     = mega ? temp[ii].objs[dind] : mw.td
        tempcorr = (exp(hnuk/tarr)-1.d) / (exp(hnuk/mw.td)-1.d)
        e_cmfn   = m0 * sarr * darr * darr * tempcorr
        e_cdf    = ccdf(e_cmfn,/cdf)
        e_cfn[ii,*] = (interpol(e_cdf,e_cmfn,mscale) > 0.) < 1.
     ENDFOR
     ;; Loop through bins to get mu, sig
     ebar = fltarr(nbin)
     FOR jj=0,nbin-1 DO ebar[jj] = stddev(e_cfn[*,jj])
     
     locfn   = (cfn - ebar) > 0.
     hicfn   = (cfn + ebar) < 1.
     clim90p = interpol(mscale,locfn,0.9) ; The 90% completeness level, lo
     clim90m = interpol(mscale,hicfn,0.9) ; The 90% completeness level, hi
     clim90pm = mean([clim90p-clim90,clim90-clim90m])
  ENDIF
  
  
  ;;===============================================================
  ;; If we wants the plots, then we makes the plots...
  IF plot THEN BEGIN
     
     pfn = './masses/plots/mass_completeness_'
     CASE ml + 2*mega OF
        1: BEGIN
           pfn += 'ml.eps'
           titstr = 'ML'
        END
        2: BEGIN
           pfn += 'mega_td'+string(sampletd,format="(I1)")+'.eps'
           titstr = 'MEGA'
        END
        ELSE: message,'Man, my liquor cabinet is getting low...'
     ENDCASE
     myps,pfn
     cdfcol = 'blk5'
     sun = cgSymbol('sun')
     
     ;; Generate differential 'completeness mass function' (DCMF)
     plothist,alog10(cmfn),xarr,yarr,bin=binw,/noplot
     yarr = [0,yarr,0] / float(ndloop)
     dx   = mean(xarr[1:*]-xarr[0,*])
     xarr = [xarr[0]-dx,xarr,xarr[n_elements(xarr)-1]+dx]
     
     ;; Plot the DCMF first
     cgPlot,10.^(xarr),yarr,psym=10,symcolor='royal blue',thick=3,charsize=1.0,$
            xr=10.^[0.25,3.05],/xst,/xlog,xmargin=[10,10],$
            yr=[0.1,max(yarr)],yst=8,ylog=1,ytickformat='exponent10',$
            xtit='Mass [M'+sun+']',ytit='N per M bin',$
            tit=titstr+' Mass Completeness Analysis'
     
     ;; Do the CDF next
     cgAxis,yaxis=1,yr=[0,1],ytit='Cumulative Distribution',/save,ylog=0,$
            charsize=1.0
     vline,/h,0.9,thick=3,color=cdfcol,linestyle=3,/log
     cgOplot,mscale,cfn,color=cdfcol,thick=10
     vline,clim90,thick=3,color=cdfcol
     cgText,clim90,0.15,charsize=0.7,$
            '  90% Completeness @ '+string(clim90,format="(I0)")+' M'+sun
     
     IF ebars THEN BEGIN
        oploterror,mscale,cfn,ebar,errcolor='red',psym=3,/nohat,errthick=3
        al_legend,/top,/right,box=0,charsize=0.7,$
                  [string(clim90,format="(I0)")+' '+cgSymbol('+-')+' '+$
                   string(clim90pm,format="(I0)")+' M'+sun]
     ENDIF
     ;; Say something about temperatures used if MEGA
     IF mega THEN BEGIN
        tdstr = sampletd ? 'Log-Normal' : $
                string(mw.td,format = "(I0,' K')")
        al_legend,/top,/left,charsize=0.8,box=0,['T!dd!n = '+tdstr]
     ENDIF
     
     ;; Re-do x-axes to make pretty
     cgAxis,xaxis=0,/xlog,/xst,xtickformat='blank_axis'
     cgAxis,xaxis=1,/xlog,/xst,xtickformat='blank_axis'
     
     
     myps,/done
  ENDIF
  
  RETURN,cfn
END
